/*
 *	Copyright (C) 2008 CRIMERE
 *	Copyright (C) 2008 Jean-Marc Mercier
 *	
 *	This file is part of OTS (optimally transported schemes), an open-source library
 *	dedicated to scientific computing. http://code.google.com/p/optimally-transported-schemes/
 *
 *	CRIMERE makes no representations about the suitability of this
 *	software for any purpose. It is provided "as is" without express or
 *	implied warranty.
 *
 *  This program is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *   You should have received a copy of the GNU General Public License
 *   along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 */
#include <src/algorithm/Evolution_Algorithm/CFL_traits.h>
#include <src/algorithm/Evolution_Algorithm/Predictor_traits.h>
#include <src/math/BoundaryCondition/ClassicalGraphDefinition.h>


#if !defined(_Particle_density_)
#define _Particle_density_

/*! \brief perform the average of the data

 given \f$\{S_i\}_i={0,..,N-1}\f$ a set of particles,
 this operator performs the transformation
 \f$S_{i+1/2} = \frac{S_i+S_{i+1}}{2}\f$, with \f$ i=0,N-1 \f$
 The case \f$i=N-1\f$ is handled by boundary conditions.
*/
template <class data>
class particle_average 
{
	OS_STATIC_CHECK(data::Dim == 1); // we may think about higher dimensions !
public :
//! \f$S_{i+1/2} = \frac{S_i+S_{i+1}}{2}\f$, with \f$ i=0,N-2 \f$.
	data operator()(const data& field) const 
	{
		OS_size nb_particules	= field.size()-1;
		return (field + (field << 1))/2.;
	}
	void update(){};
};

/*! \brief perform the inverse average of the data

 given \f$\{S_i\}_i={0,..,N-1}\f$ a set of particles,
 this operator performs the transformation
 \f$S_{i+1/2} = 2* S_{i-1}-S_{i-1/2}\f$, with \f$ i=0,N-1 \f$
 \f$ S_{-1/2} \f$ is handled by boundary conditions (determined in template data).
*/
template <class data>
class inv_particle_average 
{
public :
	OS_STATIC_CHECK(data::Dim == 1); // we may think about higher dimensions !
//! \f$ S_{i+1/2} = 2* S_{i-1}-S_{i-1/2} \f$, with \f$ i=0,N-1 \f$.
	data operator()(const data& field) const 
	{
		OS_size nb_particules	= field.size();
		if (nb_particules < 1)
			return data();
		data return_(nb_particules); 
		for (OS_int i= 0; i<nb_particules - 1;++i) 
				return_[i+1] = 2.*field[i] - return_[i] ;
		return return_;
	}
	void update(){};
};


/*!
Unclear : to rewrite completely. \f$\RR\f$

Let \f$\{ u_i \in \RR, S_i \in \RR^N \}_{i=1,..,M}\f$ be an optimal quantizer of a 
random variable u, that is 
for short a Voronoi tesselatation \f$\{\Omega_i \in \Omega \}_{i=1,..,N}\f$
determined by its location \f$\{\S_i \in \Omega \}_{i=1,..,N}\f$ such that

-
- we have \f$ {\EE}_u(\Omega_i) = u_i \f$.
- and \f$ {\EE}_u(\Omega) = \sum_{i=1,..,N} u_i = 1\f$.


This function aims to
return the discrete related density probability function of the quantizer :
we define this object as a list of
densities \f$\{v_i \in \RR, R_i \in \RR^N \}_{i=1,..,M}\f$. Let \f$\{\Delta_i\}_{i=1,..,M}\f$ 
the Voronoi tesselation induced by 
\f$\{ R_i \in \RR^N \}_{i=1,..,M}\f$.
The properties of this density probability function must be :

-
- first \f$ \int_{\Delta_i} v_i = u_i \f$.
- must be invertible.

At present time, we only provide the equi probable one dimensional case :

let \f$\{\frac{1}{N},S_i\}_{i=0,..,N-1}\f$ be a ordered one dimensional quantizer (\f$S_i < S_{i+1}\f$). 
Let us use the convention \f$S_{-1} = -\infty, S_{N} = + \infty\f$.

This quantizer defines the Voronoi tesselation
\f$\Omega_i = [S_{i-1/2}, S_{i+1/2}]\f$, i = 0,..,N

-
- with \f$S_{i+1/2} = \frac{S_i + S_{i+1}}{2}\f$ 
- and the convention \f$S_{-1/2} = -\infty, S_{N+1/2} = +\infty\f$

In this case, we define the related probability density as
\f$\{u_{i+1/2}, S_{i+1/2}\}_{i=0,..,N-1}\f$ with
\f$u_{i+1/2} = \frac{1}{ (\delta_y S)_i} = \frac{h}{S_{i+1}-S_{i}}\f$, where \f$h:= \frac{1}{N}\f$.
This object defines \f$\{\Delta_{i+1/2}\}_{i=-1,..,N}\f$  

-
- where \f$\Delta_{i+1/2} = [S_i, S_{i+1}], i= 0,..,N-1\f$
- and \f$\Delta_{-1} = [\-infty, S_0]\f$.
- and \f$\Delta_{N} = [S_N, +\-infty]\f$.

Let us verify the conditions 1) + 2)

-# \f$\int_{\Delta_i} v_i = \int_{[S_i, S_{i+1}]} \frac{h}{S_{i+1}-S_{i}} = \frac{1}{N}\f$. 
ok with the convention \f$\int_{\Delta_{-1}}\f$
-# ok (cf the inverse function below).
*/
template <class data,class result = data>
class particle_density
{
public :
	void rebuild(){
		OS_DYNAMIC_CHECK(false,"to implement");
	};
//! \f$ T_{i} = \frac{h}{S_{i+1}-S_i} \f$, with \f$ i=0,N-1 \f$.
	result operator()(const data& field) const 
	{
		OS_STATIC_CHECK(data::Dim == 1); // "I still don't know how to define a density from a particle<OS_double> field in the general case";
		if (field.size() < 1)
			return result();
//		return (1./ (double ) field.size() ) /( ( field << 1 ) - field ) + (1./ (double ) field.size() ) /(field - (field >>1) );
		return 1./(   (field<<1) - (field>>1) ) / (.5*field.size());
	};
	void update(){
		OS_DYNAMIC_CHECK(false,"to implement");
	};
};

/* !
 This operator performs the inverse transformation given by the operator
 particle density :
 We recall that the operator particle density computes
\f[
  u_{i+1/2} = \frac{h}{S_{i+1}-S_{i}}
\f]
where  \f$h:= \frac{1}{N}\f$. We are given \f$\{u_{i+1/2}\}_{i=0,..,N-1}\f$, we would like
 to compute \f$\{S_i\}_{i=0,..,N-1}\f$ verifying the previous equation.

We compute formally for \f$0 <i \le N\f$
\f[
S_{i+1} = S_i +\frac{h}{u_{i+1/2}}
\f]

We need thus an extra condition to inverse this operator. We assume
at present time \f$\sum_{i=1,..,N} S_i = 0\f$ (most of the particles
method are first moment conservative) , we will see later if a
better assumption is possible. With this assumption, we compte
recursively
\f[
    S_{i} = S_0 + \sum_{j=0,..,i-1} \frac{h}{u_{j+1/2}}.
\f]
\f[
S_0+\sum_{i=0,..,N-1} S_{i+1}  =  N S_0 + \sum_{i=0,..,N-1}
(N-i)\frac{h}{u_{i+1/2}} = 0.
\f] Thus the condition is
\f[ S_0 = -
    \sum_{i=0,..,N-1} \frac{N-i}{N^2 u_{i+1/2}}
\f]
*/
template <class data>
class inv_particle_density
{
public :
//! \f$ T_{i+1} = T-{i}+\frac{h}{S_i} \f$, with \f$ i=0,N-1 \f$.
	data operator()(const data& field) const 
	{
		OS_STATIC_CHECK(data::Dim == 1); //"I still don't know how to define the inverse of a density from a particle<OS_double> field in the general case");
		OS_size nb_particules	= field.size();
		if (nb_particules < 1)
			return data();
		data return_(nb_particules);
		return_[0] = 0;
		for (OS_int i=0; i<nb_particules;++i) {
				return_[0]  -= (nb_particules-i) /  field[i] ; 
		};
		return_[0] /= nb_particules*nb_particules;
		for (OS_int i=0; i<nb_particules-1;++i)  {
				return_[i+1]  = return_[i] + 1. /  field[i] / nb_particules; // out of range must be handled by the boundary conditions (accessor of class particle<OS_double>)
		}
		return return_;
	}
	void update(){
		OS_DYNAMIC_CHECK(false,"to implement");
	};
};



#endif